sampleSize=1000;
outputSize=1000;
lagMax=10; #Max lag for autocor
ro<-seq(length=lagMax+1, 0, 0); #autocor
Z<-seq(length=outputSize, 0, 0); #output standard gaussian
W<-seq(length=outputSize, 0, 0); #output

#generate observed samples
#X=rexp(sampleSize, rate=1);
#X=runif(sampleSize, min=0, max=1);
X=trace1;
sampleSize=length(X);

#add two extreme value
X[1]=-100000000000;
X[length(X)]=100000000000;

#sort X to generate more depedency
#X=sort(X);

#map to standard gaussian
F=ecdf(X);
Fx=F(X);
Y<-seq(length=sampleSize-2, 0, 0);
for(i in seq(2, length(X)-1, 1))
{
	Y[i-1]=qnorm(Fx[i], mean=0, sd=1);	
}
#remove first and last element of X
X=X[-length(X)];
X=X[-1];

#calculate autocor of the mapped standard gaussian
Ro=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);

#build covariance matrix (inversely), M means N-1
ro[1]=Ro$acf[1];
for(i in seq(2, lagMax+1, 1))
{
	ro[i]=Ro$acf[lagMax+3-i];
}

Ro_N_N=seq(length=(lagMax+1)*(lagMax+1), 0, 0);
attr(Ro_N_N, "dim")=c(lagMax+1, lagMax+1);
for(i in seq(1,lagMax+1,1))
{
	for(j in seq(1, lagMax+1, 1))
	{
		Ro_N_N[i, j]=ro[abs(j-i)+1];
	}
}
Ro_1_M=Ro_N_N[c(1:1),c(1:lagMax+1)];
attr(Ro_1_M, "dim")=c(1,lagMax);
Ro_M_1=Ro_N_N[c(1:lagMax+1), c(1:1)];
attr(Ro_M_1, "dim")=c(lagMax,1);
Ro_1_1=Ro_N_N[c(1:1), c(1:1)];
Ro_M_M=Ro_N_N[c(2:(lagMax+1)), c(2:(lagMax+1))];

#initialize output
mu_matrix=seq(length=(lagMax),0,0);
for(i in seq(1, lagMax, 1))
{
	Z[i]=Y[i];
	mu_matrix[i]=Y[i];
}
attr(mu_matrix, "dim")=c(lagMax,1);

temp1=Ro_1_M %*% (solve(Ro_M_M));
temp2=Ro_1_M %*% (solve(Ro_M_M)) %*% Ro_M_1;

#calculate new mu and rio
for(i in seq(lagMax+1, outputSize, 1))
{
	new_mu=temp1 %*% mu_matrix;
	new_ro=1-temp2;
	Z[i]=rnorm(1, mean=new_mu, sd=new_ro);
	
	for(j in seq(1, lagMax, 1))
	{
		mu_matrix[j,1]=Z[i-lagMax+j];
	}
}

#plot cdf
plot(ecdf(Y), col="red");
plot(ecdf(Z), col="blue", add=TRUE);

#calculate autocor
Y_acf=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);
Z_acf=acf(Z, lag.max=lagMax, type="correlation", plot=FALSE);
